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l.i  Introduction 

The  SAI/NORDA  sigma  coordinate  ocean  forecasting  com¬ 
puter  code  (sigma  code)  is  a  complex  tool  for  predicting 
ocean  temperature,  salinity,  density,  sound  speed,  and  the 
three  components  of  current  velocity  as  functions  of 
location  and  depth.  The  non-uniform  sigma  vertical  repre¬ 
sentation  is  used  to  simplify  the  treatment  of  the  severe 
topography  encountered  in  the  ocean  and  to  concentrate 
computation  levels  near  the  surface.  A  stable,  temporally 
implicit  finite-difference  representation  of  the  model 
differential  equations  allows  the  use  of  large  time  steps. 

The  basic  code  has  been  developed  under  past 
contracts  (see  references).  The  testing  of  the  model, 
however,  has  been  limited  by  the  simple  no-flux  lateral 
boundary  conditions  installed  in  the  model  during  its 
development. 

This  report  will  discuss  the  development  and  partial 
implementation  of  open  boundary  conditions  into  the  model. 

Also  presented  are  the  results  of  a  series  of  smaller 
tasks  which  have  been  undertaken  to  increase  the  code's 
usefulness.  They  include  new  mechanisms  for  data 
initialization,  new  output  options,  and  the  conversion  of 
the  code  to  the  VAX  super  minicomputer. 

Finally,  initial  results  will  be  presented  for  a 
study  of  ocean  behavior  in  the  semi-closed  basin  of  the 
Eastern  Mediterranean. 


t- 


2.0  Open  Boundary  Conditions 

The  original  form  of  the  code  used  a  very  simple 
lateral  boundary  condition;  no  flux  was  allowed.  This 
rigid  wall  approximation,  while  useful  for  testing,  is  not 
valid  for  any  ocean  region  of  interest.  All  ocean  areas 
have  significant  inflow  and  outflow. 

Two  approaches  were  taken  for  providing  open 
boundaries.  The  first  is  a  simple  condition,  which  can  be 
of  use  in  regions  where  the  flow  in  and  out  is  of  limited 
spatial  extent.  Two  examples  are  the  Gulf  of  Mexico  and 
the  Mediterranean.  This  scheme  will  be  discussed  in 
section  2.1. 

A  more  realistic  approach  requires  full  information 
about  temperature,  salinity,  and  velocities  at  each  point 
of  the  simulation  region  boundary.  To  provide  such 
information,  a  nested  run  approach  was  undertaken.  In 
this  scheme,  a  coarse  run  is  performed  over  a  large 
region.  This  run  stores  initial  and  boundary  data,  which 
is  used  by  a  second  run  whose  simulation  region  is  nested 
within  the  region  of  the  large  run.  Because  of  the  very 
complicated  spatial  representation  used  within  the  Sigma 
Code,  this  second  approach  requires  very  complex  coding. 
The  partial  implementation  of  this  second  approach  will 
be  described  in  section  2, 
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2.1  Flax  Specified  Boundaries 

In  this  approach,  a  single  flux  is  specified  for  each 
of  the  four  edges,  of  the  system.  A  function  of  lateral 
area  is  used  to  determine  how  much  of  the  total  flux 
will  be  deposited  at  each  spatial  location.  For  example, 
if  the  west  edge  of  a  simulation  region  is  all  land, 
except  for  a  narrow  strait,  virtually  all  the  flux  will 
occur  at  the  strait.  This  follows  from  the  fact  that  the 
strait  contains  all  the  lateral  area  for  that  edge  of  the 
system  (see  figure  1).  It  should  be  remembered  that  in 
the  sigma  code  there  are  a  full  set  of  computational 
levels,  even  over  land  areas.  This  scheme  ensures  that 
over  such  land  areas,  no  flux  will  be  allocated. 

The  flux  at  each  lateral  point  is  given  by 

un  «  A*f(t)*U(zk,  zb,  DINFL,  DINFLS) 

where  A  is  determined  below,  f(t)«l  -exp(-ISTP/30)  with 
ISTP  the  time  step  number,  and  D  is  the  functional 
dependence  of  the  vertical  profile.  We  define 


where 


[vi+W 


vz 

dbotbl 


DINFLS- Z _ 

(DINFLS+Z )  (1+  Z/DINFL) 


DINFL,  DBOTBL,  and  DINFLS  are  given  as  input  parameters,  z 
is  the  depth,  and  z^  is  the  bottom  depth.  DBOTBL,  through 
W,  controls  a  function  which  linearly  drops  from  z^/DBOTBL 
at  the  surface,  to  0  at  the  bottom. 

The  function  V  determines  the  actual  shape  of  the 
vertical  profile.  DINFL  sets  the  approximate  depth  of  the 
inflow.  Since  it  is  sometimes  desirable  to  have  both 
inflow  and  outflow  occuring  at  the  same  location,  V 
changes  sign  at  depth  DINFLS.  By  setting  DINFLS  very 
large  (1.E10),  no  sign  change  occurs,  and  V  is  always 
positive.  When  DINFLS  is  0.,  the  sign  change  occurs  at 
the  surface,  and  V  is  always  negative. 

The  constant  A  is  adjusted  according  to 

FLUXIN  *  -AEUan  , 

so  that  the  total  flux  at  each  edge  totals  up  to  the  input 
parameter  FLUXIN.  U  is  as  defined  above,  is  a  geometry 
factor  which  depends  on  the  edge  being  specified,  and  the 
sum  runs  over  all  points  on  the  edge. 

There  are  two  modes  of  operation  for  this  boundary 
condition.  The  first  is  unidirectional  flow.  The  user 
sets  DINFLS  large  O1.E10),  and  sets  DBOTBL,  DINFL,  and 
FLUXIN  for  the  desired  flow. 

The  second  mode  is  bidirectional  flow.  In  this  case, 
the  user  can  set  the  ratio  of  inflow  to  outflow,  DINFLR. 
The  code,  by  an  iterative  scheme,  adjusts  one  of  the 


other  variables  to  obtain  the  correct  ratio.  If  DINFLR>0, 
DINFLS  is  adjusted.  DINFLR<0  is  a  way  to  adjust  DINFL 
instead.  For  best  results,  a  good  guess  for  DINFL  and 
DINFLS  should  be  made.  By  running  the  code  with  NSTEP«0, 
the  final  value  of  DINFLR  can  be  compared  to  that 
requested. 

The  input  variables  are: 

FLUXIN(i)  total  flux  (outward)  for  this  edge, 

DBOTBL(i)  scaling  factor  for  linear  profile, 

DINFL(i)  depth  of  inflow, 

DINFLS (i)  depth  of  sign  change,  and 

DINFLR (i)  ratio  of  inflow  to  outflow. 

Here  i  determines  which  edge  of  the  system  according  to 

i*l  west  edge 

2  east  edge 

3  south  edge 

4  north  edge 

The  adjustment  of  DINFL  or  DINFLS,  and  the  determination 
of  A  is  performed  in  new  subroutine  FLUXST.  The  flux 
calculated  by  this  scheme  is  added  to  the  velocity  arrays 
in  DVBND. 


2.2  Open  Boundary  Conditions  Based  on  a  Nested  Run 

A  nested  run  consists  of  a  fine  resolution  run  whose 
simulation  region  is  located  within  the  simulation  region 
of  a  coarse^  larger  area  run.  (See  Figure  2).  There  are 
four  aspects  of  such  a  simulation.  One,  a  coarse 
simulation  is  performed.  This  run  writes  out  both  data 
for  initialization  of  the  fine  run,  plus  boundary  data  on 
the  embedded  boundaries.  Second  the  fine  run  performs  a 
complex  interpolation  from  this  data  to  obtain  an  initial 
condition.  Third#  the  fine  run  performs  a  complex 
interpolation#  at  each  time  step#  to  find  appropriate 
values  along  its  boundary.  Fourth#  these  boundary  values 
are  used  in  a  suitable  boundary  condition.  Each  of  these 
steps  will  be  described  below. 

2.2.1  Coarse  Run 

The  initialization  data  written  out  by  the  coarse  run 
consists  of  temperature#  salinity#  pressure  and  velocities 
over  a  region  which  encompasses  the  fine  region.  All  the 
coarse  data  is  not  needed#  and  it  is  not  written  out. 
There  may#  however#  be  more  than  one  fine  region  nested 
within  a  particular  coarse  region.  The  data  for  each  fine 
region  is  written  to  a  different  file. 

The  number  of  such  subregions  is  specified  by  input 
parameter  NREGON.  For  each  of  these  NREGON  subregions# 
the  location  of  the  boundaries  is  specified  by: 
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SAVLLO(i)  west  edge  of  region  i  in  degrees 
SAVLHI(i)  east  edge  of  region  i  in  degrees 
SAVPLO(i)  south  edge  of  region  i  in  degrees 
SAVPHI(i)  north  edge  of  region  i  in  degrees. 

The  coarse  run  uses  these  input  variables  in 

subroutine  REGNST  to  establish  the  areas  on  its  grid  that 

will  have  to  be  written  out  for  each  of  the  subregions. 

It  also  determines  the  locations  of  its  grid  lines  which 

lie  on  either  side  of  the  embedded  boundaries,  for  both  T- 

S  and  U-V  data,  again  for  each  subregion. 

At  each  time  step,  subroutine  SAVREG  is  called  to 

write  out  the  appropriate  data.  The  writing  of  data  is 

controlled  by: 

SAVSTR(i)  time  step  of  initial  write  for  region  i. 
SAVlNC(i)  modulus  on  time  step  to  perform  write  for 
region  i. 

On  only  the  initial  write,  data  for  initialization  is 
written  out.  On  this,  and  all  subsequent  writes,  boundary 
data  is  written.  Each  subregion  has  its  own  file. 
Therefore,  the  location  and  frequency  of  the  data  may  be 
completely  different  for  each  subregion.  The  file  unit 
number  is  70  +  i,  where  i  is  the  region  number. 

Subroutine  RWBND  is  used  to  read  or  write  the 
boundary  data. 

2.2.2  Initialisation  of  Final  Run 

The  data  written  out  by  the  coarse  run  must  be 
interpolated  onto  the  fine  mesh  of  the  nested  run.  This 
is  complicated  by  the  complex  spatial  representation  of 


the  sigma  code.  As  can  be  seen  in  Figure  2,  the  lateral 
mesh  spacing  may  be  nonuniform  and  arbitrary.  As  shown  in 
Figure  3,  the  vertical  spacing  is  nonuniform  and  depends 
both  on  the  topography,  and  the  number  of  vertical  levels. 
In  a  typical  simulation,  no  two  grid  nodes  are  located  at 
the  same  depth. 

The  initialization  in  the  sigma  code  is  performed  by 
one  of  the  "CASE”  routines.  For  the  nested  run  the  input 
parameter  CASE  is  set  to  8.  This  causes  subroutine  CASES 
to  be  called  to  perform  an  initialization  from  coarse  run 
data.  Storage  must  be  carefully  handled  so  that  no  excess 
array  space  is  allocated.  To  do  this,  CASE8  was  designed 
as  a  dummy  routine  to  handle  storage  for  subroutine  CASE8A 
where  the  actual  interpolation  is  performed. 

In  CASE8A  an  interpolation  is  performed  using  the  8 
surrounding  data  points  from  the  coarse  data.  Weight 
factors  are  determined,  and  then  the  new  value  is 


calculated  from. 


I  «(i)xc(i) 

xf =  - 

i  u(i) 
i=l 


Because  of  the  representation,  the  eight  weights  have 
different  values  for  each  of  the  fine  mesh  grid  points. 


Additional  problems  arise  because  of  the  use  of  real 
topography.  In  Figure  3,  point  A  is  a  possible  location 
of  a  fine  mesh  point.  It  can  use  the  coarse  value  to  the 


left,  but  the  coarse  value  to  the  right  does  not  lie  deep 
enough  for  point  A,  and  it  cannot  be  used.  The 
interpolation  scheme  is  therefore  modified  to  exclude  the 
point  to  the  right,  and  the  weights  for  the  points 
actually  used  are  corrected. 

Point  B  is  the  location  of  a  second  type  of  problem. 
At  this  location  no  coarse  data  is  available  at  the  depth 
of  point  B.  For  this  problem  a  different  procedure  is 
used.  Subroutine  CASE8B  is  called  to  perform  a  spline 
extrapolation  of  the  surrounding  data  to  obtain  values  at 
the  correct  depth.  Then  a  four  point  lateral 
interpolation  of  these  extrapolated  values  is  used  to 
obtain  a  reasonable  result.  This  involved  procedure  is 
only  used  for  temperature  and  salinity  values.  For  a 
velocity  point  at  location  B,  when  no  coarse  data  are 
available,  the  velocity  is  set  to  zero.  This  ensures  that 
no  artifical  flows  will  result. 

It  should  be  noted  that  different  weights  are 
required  for  temperature-salinity  (T-S)  and  velocity  data, 
since  T-S  data  locations  are  spatial  shifted  from  velocity 
locations. 

2.2.3  Determination  of  Boundary  Value  Data  for  Fine  Run 

At  each  step  during  the  fine  run,  the  values  of 
temperature,  salinity,  the  velocity,  and  the  pressure  must 
be  determined  at  the  location  of  the  physical  boundaries 
of  the  smaller  system.  These  values  are  calculated,  and 


put  into  a  set  of  two-dimensional  arrays.  These  arrays 
are  used  by  the  boundary  conditions  in  determining  the 
values  of  the  relevant  quantities  at  the  guard  cells.  It 
should  be  noted  that  these  boundary  values  are  not  used 
directly,  but  are  used  in  the  application  of  the  boundary 
conditions. 

Data  is  written  out  on  both  sides  of  the  boundary  by 
the  coarse  run.  Since  T-S  and  U-V  data  are  spatial 
displaced,  these  surrounding  boundary  planes  are  located 
differently  for  the  two  types  of  data. 

Subroutine  EXBND  is  used  to  determine  the  values  of 
data  on  the  boundaries.  As  in  the  use  of  CASE8,  EXBND 
acts  as  a  storage  manager  for  EXBNDA,  where  the  actual 
interpolation  is  performed.  The  interpolation  scheme  is 
the  same  as  that  used  in  CASE8A,  with  subroutine  CASE8B 
used  here  also  for  points  out  of  range. 

2.2.4  Implementation  of  Rested  Boundary  Conditions 

Our  planned  implementation  of  nested  boundary 
conditions  is  based  on  a  major  extension  of  the  analysis 
in  §  4.4  of  the  prior  study  (Roberts,  1982).  In  that 
section  we  studied  the  reflection  and  transmission 
properties  of  different  sets  of  open  boundary  conditions 
applied  to  the  one-dimensional  wave  equation,  and 
demonstrated  that  the  spurious  waves  could  be  reduced  to 


a  minimum  by  imposing 


cun  +  p 

on  the  boundary,  where  p  is  the  scaled  pressure,  un  is  the 
outward  flow,  and  c  is  the  local  phase  speed  of  the  wave. 
The  success  of  this  method  follows  from  the  fact  that  the 
outward  and  inward  propagating  wave  modes  have 

cun  *  ±p, 

respectively. 

In  the  prior  report  we  reached  the  conclusion  that  a 
generalization  of  this  boundary  condition  to  the  three- 
dimensional  sigma  code,  with  rotation,  diffusion, 
stratification,  and  bottom  topography  effects,  in  addition 
to  the  surface  gravity  waves,  was  impractical.  In  §  5  we 
outlined  a  simpler  set  of  boundary  conditions  based  on  the 
normal  flow  and  on  conservation  requirements. 

However,  we  now  believe  that  for  most  situations  of 
interest,  spurious  surface  and  internal  gravity  waves 
generated  at  the  open  boundaries  by  poorly  chosen  boundary 
conditions  will  pose  the  greatest  problem  in  nested 
computations.  We  therefore  plan  to  implement  a 
generalization  of  the  above  boundary  condition  to  the  full 
three-dimensional  sigma-coordinate  system  of  equations. 

Our  boundary  conditions  are 


r 


Qn  =  und  +  F(i  V  +  eg(z»zb) 


6nu||  +  a(u||  -  u||d>  =  0 


<SnTn  +  a(Tnn-  Td)  =  0 


6nSn  +  (Snn-  Sd)  =  0 

Here 

un  and  Uj,  denote  the  velocity  components  outward  and 
parallel  to  the  boundary, 

fn  and  f^ denote  averages  of  adjacent  mesh  values  in  the 
same  directions,  for  any  variable  f, 

f (3  denotes  data  values  obtained  by  interpolation  from  the 
previous  coarse  run,  at  the  appropriate  positions  and 
times, 

a  is  a  positive  function  of  un6x/KH,  and  is  zero  for  large 
values  (passive  outflow  condition)  and  very  large  for 
negative  values  (value  imposed  on  inflow), 

6nTn  is  the  difference  over  two  mesh  intervals  (outside 
minus  inside)  divided  by  two, 

¥nn  is  the  mean  of  the  outside  and  inside  values,  ignoring 
the  boundary  value  (a  non-standard  notation), 
g  is  (z-zb)/(l+z/zi) 

3  is  zero  or  is  chosen  so  that  the  total  normal  is  an 
imposed  value  obtained  from  data  on  the  movement  of  the 
free  surface. 

F  is  an  operator  designed  to  separate  the  internal  wave 
modes,  and  divide  each  by  an  appropriate  wave  speed  c;  the 
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simplest  operator  option  is  to  divide  by  a  single  imposed 
c  value. 

These  boundary  conditions  replace  our  previous 
conditions 

0n.0 

6nu||  -  0 

Snin  *  ° 

6  Sn  =  0 
n 

For  our  Mediterranean  simulation  we  imposed  a  nonzero  U 
distribution;  u  has  otherwise  been  zero. 

2.2.5  Status  of  implementations 

The  first  three  steps  have  been  implemented/  and 
preliminary  tests  performed.  Coarse  runs  have  been  used 
to  write  out  initial  and  boundary  data.  Fine  runs  have 
read  in  initial  data,  and  have  interpolated  this  data  to 
get  an  initial  condition  for  the  system.  In  addition,  the 
boundary  data  has  been  interpolated  to  obtain  values  of 
the  relevant  quantities  on  the  boundaries. 

There  has  been  insufficient  time,  however,  to 
implement  the  algorithms  of  Section  2.2.4.  These 
algorithms  will  go  into  subroutines  TSBND  and  UVBND,  and 
will  use  the  boundary  value  data  which  now  exists  in  the 


Subroutines  which  have  been  changed  and  new 


subroutines  are  listed  in  Appendix  A. 


3.0  External  Data  for  Initialisation 


The  quality  of  a  simulation  is  very  dependent  on  the 
quality  of  the  data  used  for  the  initialization  of  the 
run.  It  is,  therefore,  very  important  to  use  the  best 
quality  data  available  for  initialization.  OTS  and  EOTS 
data  is  of  better  quality  than  climatology,  but  such  data 
does  not  go  deep  enough  to  initialize  the  entire 
simulation  region.  The  best  solution,  therefore,  is  to 
smoothly  merge  two  data  sets  to  provide  OTS  or  EOTS  where 
available,  and  climatology  where  not.  The  transition 
between  the  two  data  sets  must  be  smooth,  or  numerical 
noise  will  result. 

If  Te  it,  OTS  or  EOTS  data,  and  Tc  is  climatology 
data,  the  procedure  is  as  follows.  For  depths  above  the 
lowest  Tfi  point  (z«d) 

T  -  Te(2) ,  (zid) 

where  Tg  at  arbitray  values  of  z  are  determined  from  a 
spline  fit.  For  z>d, 

T  -  Tc(z)  +  ATf (z)  +  AT'g(z)  . 

Here  Tc(z)  is  again  obtained  from  a  spline  fit,  and 
AT  «  TE  -  Tc  at  z  «  d 

AT'  •  t»e  -  t'c  at  z  =  d 

The  functional  forms  of  f  and  g  are  given  as 

f(z)  »  |exp  (l“|)r 
g(z)  ■  (z-d)  exp(l-f). 
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At  z  ■  d  these  functions  result  in 

f (d)=l  f  f'(d)«0, 
g(d)-0  ,  g« <d)«l. 

The  implementation  of  this  capability  involved  two  sets  of 
changes.  First  the  data  tape  handling  program  TAPER2  had 
to  be  slightly  modified  to  create  a  data  file  with  both 
OTS  (or  EOT S)  and  climatology.  Secondly,  two  new  CASE 
routines,  CASE6  and  CASE7  were  written  to  read  in  the 
combined  data,  and  to  perform  the  merging.  CASE6  merges 
OTS  and  climatology,  while  CASE7  handles  EOTS  and 
climatology. 

In  practice,  two  TAPER  runs  are  used  to  produce  the 
combined  input  data  set.  The  listing  in  Appendix  B  is  an 
example  which  shows  the  procedure.  A  TAPER2  run  is  used 
to  read  in  the  appropriate  OTS  (or  EOTS)  data.  Then  a 
TAPER3  run  is  used  for  the  climatology  data.  The  output 
of  TAPER3  is  appended  onto  the  output  file  from  TAPER2  to 
form  the  required  input  data  file. 

In  the  sigma  code  run,  one  simply  assigns  the  data 
file  produced  by  the  TAPER  runs.  Then  the  CASE  input 
parameter  is  set  to  6  for  OTS,  or  7  for  EOTS,  and  the 
proper  initialization  is  performed. 


4.0  Enhanced  Output  Options 


c 


© 


c 


It 


The  sigma  code  has  an  extensive  set  of  output 
options.  Graphics  output  can  be  produced  of  all  relevant 
data  fields,  on  either  plotting  devices,  or  as  printer 
plots.  Several  desired  enhancements  in  the  existing  sigma 
code  output  were  identified,  however,  and  they  are  the 
subject  of  the  current  section. 

The  geostrophic  velocity  is  the  velocity  that  results 
from  the  balance  between  the  horizontal  pressure  gradient 
and  the  horizontal  component  of  the  coriolis  acceleration. 
It  provides  useful  information  about  the  state  of  the 
physical  system  being  simulated.  Therefore,  plots  of  the 
geostrophic  velocity  were  added  as  an  output  option. 

The  geostrophic  velocity  is  calculated  during  the 
vertical  sweep  through  the  mesh  by  DVDOTB.  Rather  than 
allocate  an  additional  three  dimension  array  to  hold  this 
information,  each  two  dimensional  slice  is  written  to  a 
disk  file  by  UVDOTB  as  calculated.  At  the  end  of  the 
sweep,  new  subroutine  OVGEO  reads  this  information  into  a 
three  dimensional  scratch  array,  from  which  it  is  written 
onto  the  plot  file. 

Subroutine  OCPLOT,  in  the  plotting  package,  was 
modified  to  create  plots  from  this  newly  available  data. 
A  typical  output  plot  is  shown  in  Figure  4. 

In  addition  to  this  new  capability,  three 


enhancements  were  made  to  the  existing  code.  First,  on 
restarts,  because  of  the  way  time  step  numbering  was 
handled,  plot  specification  did  not  work  properly.  This 
was  corrected  by  a  complete  revision  of  time  step  handling 
on  restarts. 

Second,  the  orientation  of  the  plots  of  vertical 
profiles  was  awkward,  and  did  not  follow  the  conventional 
orientation  for  such  plots.  This  has  now  been  corrected. 

Third,  up  to  now,  certain  quantities  were  only 
available  as  printer  plots.  Now  such  information  can  be 
displayed  on  plotting  devices.  Sample  plots  of  the 
barotropic  stream  function  are  shown  in  Figure  5  and  6. 

The  required  changes  to  the  sigma  code  are  included 
in  Appendix  A.  Those  subroutines  in  the  plotting  program 
which  required  change  are  listed  in  Appendix  C. 

The  capability  to  plot  velocities  at  a  constant  depth 
is  available  by  modification  of  the  TSPOP  post  processor. 
In  the  main  program,  the  lateral  averaging  of  the  bottom 
topography  is  removed.  In  addition  the  calls  to  MESHST 
are  changed  so  that  the  U-V  positions,  rather  than  the  T-S 
positions  are  produced.  The  resulting  code  is  given  in 
Appendix  D. 


5.0  VAX  Version  of  Code 


The  usefulness  of  the  sigma  code  increases  when  it 
becomes  available  on  additional  computer  systems.  Because 
of  the  wide  spread  availability  of  the  VAX  11/780  super 
minicomputer,  it  was  decided  to  install  the  sigma  code  on 
this  machine. 

The  code  has  been  installed  on  the  VAX  and  a  sample 
output  plot  from  a  simple  test  run  in  shown  in  Figure  7. 
The  code  runs  much  in  the  same  way  as  it  does  on  the  TI- 
ASC  computer.  The  VAX  files  used  by  sigma  code  are  listed 


below: 

Unit 

Number 

Name 

Use 

Preexisting 

1 

SIGRUNLOG 

Log  File 

Yes 

2 

SIGBTOPOG 

Bottom  Topography 

Yes 

3 

S1GTSPLOT 

T-S  Plot  Data 

4 

SIGBUOYIN 

Buoyancy  Coefficients 

Yes 

5 

SIGINPUTD 

Input  Data 

Yes 

6 

SIGPRINTF 

Print  File 

10 

SIGUVPLOT 

U-V  Plot  Data 

14 

SIGBDOYOU 

Buoyancy  Coefficients (output) 

18 

SIGINITIA 

Initialization  Data 

Yes 

19 

SIGFORCEG 

Forcing  Data 

Yes 

26 

SIGQCHEKl 

Qcheckl  Print  File 

27 

SIGQCHEK2 

Qcheck2  Print  File 

40 

SIGRESTRT 

Restart  Data 

Yes 

42 

SIGDUMPFL 

Dump  File 

80 

SIGWORKSP 

Work  Space 

These  files 

noted  as  preexisting,  must 

be  available 

before  the  run,  as  needed. 

Because  it  was  not  known  what  graphics  software  would 
be  available  on  a  particular  VAX,  no  graphics  device 


output  is  currently  installed.  However,  printer  graphics 
are  installed,  and  working  as  shown  in  Figure  7. 

An  estimate  of  relative  timing  is  difficult,  as  VAX 
timing  depends  on  many  factors,  including  system  loading. 
A  crude  estimate  shows  that  the  VAX  version  of  the  code 
runs  about  a  factor  of  120  slower  than  the  ASC  version. 
This  is  about  the  ratio  that  was  expected. 


6.0  Mediterranean  Tests 

The  sigma  code  was  setup  to  perform  studies  of  the 
eastern  basin  of  the  Mediterranean  Sea.  The  process 
involved  routine  setup  tasks  as  well  as  two  code 
modifications.  The  setup  procedure,  as  well  as  the  code 
modifications,  will  be  discussed  below. 

6.1  Setup 

The  region  was  established  as  10°-37°E,  and  32° -38° 
N.  This  includes  the  eastern  boundary,  but  cuts  off 
several  features  at  the  north  and  south.  Although  the 
bottom  topography  was  available  at  10'  intervals,  a  20' 
spacing  was  chosen  to  reduce  memory  requirements  for  this 
problem.  The  numbers  of  lateral  grid  points  are: 

(37-10) *3+2=83  in  the  east-west  direction,  and 
(38-32) *3+2=20  in  the  north-south  direction. 

Five  vertical  levels  were  used.  The  three  code  parameters 
IF,  JF,  and  KF  were  reset,  and  a  recompilation  of  relevant 
routines  was  performed. 

The  topography  was  created  from  the  10'  data  file  by 
use  of  the  program  tested  below  (Note:  All  programs  and 
listings  referred  to  in  this  section  can  be  found  in 
Appendix  E).  A  rigid  wall  was  imposed  on  the  topography 
at  all  lateral  boundaries,  exept  for  the  Strait  of  Sicily, 
where  the  inflow-outflow  occurs. 


Since  the  forcing  and  initialization  data  lies  on  a 
polar-stereographic  (PS)  grid,  the  PS  limits  of  the 
simulation  region  were  determined.  PS  data  was  found  to 
be  required  in  the  grid  whose  indices  ran  from  1*45-50, 
and  J*31-40. 

An  initialization  file  was  created  from  climatology 
using  the  TAPER3  program.  This  was  adequate  for  testing 
purposes.  For  a  better  initialization,  the  procedure  of 
Section  3  of  this  report,  may  be  used. 

The  forcing  data  was  extracted  from  the  data  tape 
"SAIATM2",  by  use  of  the  TAPER  program.  As  shown  in  the 
listing,  the  data  starts  on  January  7,  1977.  A  six  hour 
interval  was  chosen. 


6.2  Code  Modifications 

The  original  formulation  for  flux  specified 
boundaries  allowed  only  inflow  or  outflow  at  each  edge. 
The  full  formulation,  as  described  in  Section  2.1,  was 
implemented  to  allow  bidirectional  flow,  as  required  in 
the  Mediterranean  Tests. 

The  second  addition  to  the  code  was  required  because 
of  a  fundamental  difficulty  in  the  code.  The  code  solves 
for  the  time  evolution  of  a  general  variable  x  by  the 
process, 


x  *  A(x) 

xn+l  .  xn  +  f ( x , dt ) . 


c 


A  is  a  representation  of  the  model  differential  equations, 
and  F  is  a  fixing  (stabilization)  operator.  For  a  simple 
explicit  scheme  F(x,dt)  *=  xdt.  However,  such  a  scheme  is 
unstable,  and  very  complicated  procedures  have  been  used 
to  develop  appropriate  formulations  for  F,  for  the 
equations  to  be  solved  here  (see  references). 

Since  the  Sigma  Code  solves  a  three  dimensional 
problem  in  a  minimum  of  computer  memory,  careful  attention 
was  paid  to  the  code  architecture.  The  problem  is  solved 
in  slices,  with  intermediate  quantities  overwritten  as 
each  level  is  solved.  This  dictates  the  implementation  of 
the  fixing  operator.  Specifically,  lateral  fixing  is 
performed  during  the  vertical  sweep,  before  the  vertical 
fixing,  which  requires  a  knowledge  of  all  the  vertical 
levels. 

The  problem  arises  because  the  lateral  fixing 
propagates  the  unstable  solution  horizontally,  before  the 
vertical  fixing  has  modified  the  result.  The  vertical 
fixing  can  subsequently  stabilize  the  solution  at  each 
point,  due  to  the  physics  at  each  point.  It  cannot, 
however,  stabilize  the  part  of  the  instability  which  has 
been  horizontally  propagated. 

In  the  past,  this  problem  has  been  circumvented  by 
turning  off  the  lateral  fixing,  thus  avoiding  the 
horizontal  propopagation  of  the  instability.  Because  of 
the  smaller  lateral  mesh  spacing  required  in  the 


Mediterranean  tests,  however,  the  lateral  fixing  must  be 
on. 

One  solution  is  to  add  a  preliminary  vertical 
diffusion  model,  before  the  horizontal  fixing.  This  can 
take  a  number  of  forms.  The  simple  solution  implemented 
here  is  to  limit  x  (i.e.  u),  before  the  horizontal  fixing 
process.  A  new  routine  DDLIM,  modeled  on  routine  ULIM  has 
been  written.  It  sets 

.  v*u 

u  =  v+Taf  , 

where  V^DOMAXN'DT^zjj.  DUMAXN  is  under  input  control,  DT 
is  the  time  step,  and  z^  is  the  local  depth. 

6.3  Results 

A  preliminary  test  run  was  performed.  This 
simulation  is  intended  to  demonstrate  that  the  code  is 
properly  setup  to  perform  studies  of  the  eastern 
Mediterranean.  Mo  attempt  was  made  to  model  actual 
physical  processes.  Rather  this  run  establishes 
procedures  for  modeling  the  eastern  Mediterranean. 
Additionally,  in  the  process  of  performing  this  test, 
problems  with  the  code,  when  applied  to  this  region,  were 
identified  and  corrected. 

It  should  be  noted  that  in  the  figures  which  follow, 
problems,  such  as  inadequate  labeling  of  contour  lines, 
are  problems  in  the  Disspla  plotting  package.  They 
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result.  In  part,  from  the  fact  that  the  ASC  computer 
center  is  using  an  out  of  date  version  of  Disspla. 

Figure  8  is  a  plot  of  depth  values  in  the  region. 
The  values  in  the  interior  are  the  values  obtained  from 
the  data  file.  The  boundaries  demonstrate  the  rigid  wall 
that  was  articif ically  imposed.  The  one  open  spot  is  at 
the  Strait  of  Sciliy.  Because  the  labeling  in  the  Disspla 
plot  is  sparse,  a  printer  plot  of  the  same  data  is  shown 
in  Figure  9.  Since  each  contour  is  labeled  by  a  letter, 
actual  values  can  be  read  off  this  plot. 

A  two  week  interval  was  modeled  using  56  time  steps 
of  6  hours  each.  A  great  number  of  output  plots  can  be 
produced  from  the  data  files  which  were  generated  during 
this  run.  Indeed,  there  is  no  limit  to  the  data 
representations  available.  Most  major  quantities  can  be 
displayed  on  vertical  slices  of  arbitrary  great  circle 
arcs.  Additionally,  horizontal  slices  can  be  displayed  at 
arbitrary  depths.  The  plotting  programs  perform  the 
required  interpolations,  so  that  the  spatial 
representation  in  the  code  does  not  impose  limits  on  how 
the  data  may  be  displayed.  A  representative  sample  of 
possible  output  plots  follows. 

The  general  shape  of  the  vertically  integrated  flux 
as  shown  in  Figure  10,  with  actual  values  appearing  in  the 
printer  plot  of  figure  11.  The  remainder  of  the  plots 
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show  quantities  on  the  great  circle  with  endpoints  at 
10°  E,  35°  N  and  37°  E,  35°  N.  The  temperature  is  displayed 
in  Figure  12.  This  figure  has  the  depth  adjusted  so  that 
the  entire  vertical  extent  of  the  simulation  is  shown. 
Figure  13  is  the  same  display  with  the  maximum  depth 
adjusted  to  600  m.  Such  a  plot  allows  closer  examination 
of  surface  dynamics.  The  same  set  of  plots  for  salinity 
are  shown  in  Figures  14  and  15.  Once  again  the  entire 
vertical  domain,  and  a  closer  look  at  the  surface  are 
shown . 

A  similar  set  of  plots  were  also  generated  from  the 
velocity  data.  Figure  16  is  a  plot  of  u  velocity.  For 
this  slice,  a  positive  u  velocity  is  a  velocity 
approximately  to  the  right.  Since  this  is  a  display  on  a 
great  circle  arc,  rather  than  an  arc  of  constant  latitude, 
the  velocity  is  not  exactly  in  the  plane  of  the  plot. 
Figure  17  is  a  plot  of  v  velocity.  Again,  for  this  slice 
which  runs  approximately  east-west,  a  positive  v  velocity 
is  a  velocity  roughly  into  the  page. 

The  last  two  plots,  Figures  18  and  19,  demonstrate 
the  new  ability  of  the  code  to  plot  the  two  components  of 
the  geostrophic  velocity.  The  comments  made  in  relation 
to  the  velocity  also  apply  here. 


6.4  Conclusion 

The  Sigma  code  is  now  working  in  the  eastern  basin  of 
the  Mediterranean.  Externally  supplied  data  has  properly 
been  integrated  into  the  code  for  both  bottom  topography 
and  initialization.  In  addition,  a  combined  inflow- 
outflow  condition  at  an  open  strait  is  working. 

Actual  simulation  studies  will  involve  longer  runs. 
A  detailed  examination  of  the  possible  output  displays, 
made  in  light  of  the  topography  of  the  region,  will 
provide  an  insight  into  the  processes  at  work.  As  such 
studies  are  performed,  it  may  become  evident  that  changes 
in  the  externally  supplied  data  are  necessary.  These  may 
include  the  addition  of  new  data  (e.g.  wind  forcing),  or 
the  use  of  better  data  than  is  currently  used  (e.g.  EOTS 
rather  than  climatology). 
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1  -  - 

LANGUAGE: 

USER  1  (FORMATION: 

1 

/  CON 

TWODATA 

1 

2 

/  CON  MOWS: 

TWOOATA 

2 

m 

3 

/  CON 

TWODATA 

3 

4 

/CON  GET  IS  LIKE  ASG 

TWOOATA 

4 

8 

5 

/CON  RGET  IS  LIKE  REL  AM)  T)€N  ASG 

TWODATA 

5 

■/  ] 

6 

/CON  PUT  IS  LIKE  CATV 

TWODATA 

6 

.-•< 

7 

/CON  SEE  IS  LIKE  FOSYS 

TWOOATA 

7 

•• 

8 

/  CON 

TWODATA 

8 

9 

/CON  PT  IS  THE  PATH  FOR  TFE  TAPER  OBJECT  LIBRARY 

TWOOATA 

9 

• 

10 

/  CON 

TWODATA 

10 

11 

/  UNIT  MIN=2 

TWODATA 

11 

f 

12 

/  RGET  OBJLI8.PT/08JLIB 

TWODATA 

12 

13 

/  REL  SYS.LNOO 

TWODATA 

13 

14 

'  UK 

TWODATA 

14 

15 

LIBRARY  QBJLIB 

TWODATA 

15 

* 

m 

16 

INCLUDE  TAPER2 

TWODATA 

16 

17 

/  REL  FT31F001 . FT32F001 >  FT33F001 . R34F001 . FT35FOO 1 . FT36F001 

TWODATA 

17 

> 

18 

/  REL  FT41F001,FT42F001.FT43F001.FT44f001.FT45FOOI.FT46F001 

TWODATA 

18 

19 

/  REL  FT51F001 . FT52F001 . FT53F001 . FT54F001 . FT55F001 >  FT56F00 1 

TWOOATA 

19 

20 

/  Ra  FT71F001.FT72F001 

TWODATA 

20 

21 

/  RGET  FT61FOOt.DOO/NAVY/NOROA/UARNAl/FNOCTAPE/SAIQCNl 

TWODATA 

21 

22 

/  SET  N=l 

TWODATA 

22 

< 

23 

/  REL  FT06F001 

TWOOATA 

23 

-  A 

24 

/  FO  FT06F001.BANB=2/20/2 

TWODATA 

24 

25 

/  FO  FT71F001.BAMW/50/4 

TWOOATA 

25 

’  ■  ■■■■! 

26 

/  FXOT  QPT=(I).LTP=<99.99.N1.AD0NEHI24K.CPTIHE=6Q00 

TWOOATA 

26 

27 

60NE 

TWODATA 

27 

'  '  ■/; 

28 

IS=05.  IE=28,  JS=13.  JE=51. 

TWOOATA 

28 

■  -i 

LA 

29 

YEAR=76.HONTH:lO.DAY=29.HOUR*Oi 

TWOOATA 

29 

30 

0T=24 . . NT1 1 . NCATST=7 . 

TWOOATA 

30 

r  . 

31 

SIGNAD=T. 

TWODATA 

31 

32 

CATN0='B10+' . 'P14+' , 'PI5+' . 'P16+' , 'PI7v . 'PI8*' . 'PI9+' . 

TWOOATA 

32 

33 

6EW 

TWOOATA 

33 

34 

/  CON 

TWOOATA 

34 

35 

/  CON  OUTPUT  ON  FT71F001 

TWOOATA 

35 

\9 

36 

/  CON 

TWODATA 

36 

»-  < 

i 

37 

/  IF  TERH.MLO.ERR 

TWODATA 

37 

38 

/  SEE  N.FT06F001,NANE=FIl£60NE 

TWOOATA 

38 

39 

/  RGET  0BJLIB.pt/08JLIB 

TWODATA 

39 

-  * 

40 

/  REL  SYS.LNOO 

TWODATA 

40 

41 

/  UK 

TWOOATA 

41 

r 

42 

LIBRARY  OBJLIB 

TWODATA 

42 

.  ■% 

!C 

43 

INCLUDE  TAPER3 

TWODATA 

43 

i 

44 

/  REL  FT31F001  >  R32F001 .  R33F001 .  FT34F001 .  FT35F001 .  FT36F001 

TWODATA 

44 

45 

/  REL  R41FOOl.n42F001.FT43F001.n44F001.FT45FOOt.R46F 001 

TWODATA 

45 

46 

/  Ra  FT51F001 .  FT52F001 .  F753F001  .FT54F001 .  FT53F001 .  FT56F001 

TWODATA 

46 

47 

/  CON 

TWOOATA 

47 

48 

/  CW  TAPER3  WRITES  ON  R72  ->  CHANGE  SO  THAT  *  WILL  APPEND 

TWODATA 

48 

it 

49 

/  cm 

TWOOATA 

49 

p 

50 

/  Ra  FT72F001 

TWODATA 

50 

51 

/  RENAME  FT71F001.FT72F001 

TWOOATA 

51 

52 

/  FO  FT72F00I.P0SNB0 

TWODATA 

52 

53 

/  RGET  FT61F00L , AFFIL/ IND/SAI /SEFTJ1 /USER/fCDN/CLIHAT 

TWOOATA 

53 

i 

54 

/  SET  N«1 

TWODATA 

54 

< 

t 

/ 

> 

f 
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55 

/  PEL  FT06F001 

TNOOATA 

55 

56 
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Main  Program  for  UVPLOT 


Win  I  MCI  AAtM  IS  fOlliO  »M  sn  FILE 
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02/15/04 


BEAia  OR  02/I5/D4  AT  It » 1*03 


LAST  UPDATED  ON  02/15/04  AT  11 « 1*03  BY  SHB  \««ION  3.27 


USED  nrOMATIONs 


FXL  LIWARWMJLH 
PNOORAH  BUTTON 


tmmtt  solution  domain  itmmtt 
nmmmmtntmmttmtntm 
mmutmttutmtmmmtum 
mmmimmmmmtmmtm 


to  t 


T(l)  T(2) 


T(IFPU 


14  t 


1(1)  1(2) 


Tt€  ACTUAL  SOLUTION  00NA1N  NAS  IF-2  INTERVALS  AN)  RUNS  FR® 

T (2)  TO  T(IF>.  (Til)  AND  T(IFPl)  ARE  EXTRA  POINTS.) 

I  IS  STAGGERED  WT  T  AND  IS  OEFDO  FROH  A  HALF  CELL 
TO  THE  LEFT  OF  THE  DOMAIN  TO  A  HALF  ICSH  CELL  ON  THE  RIGHT. 

ME  PASS  TO  ROUT  IIC  INTRP  ACTUAL  LOCATIONS  IN  DEGREES.  UtCAE 
NE  ISO  TIC  DEPTH  (Z).  THAT  IS  AN  AREA  HHIQl  IS  LARGER 
THAN  ITC  DOMAIN  BY  A  HALF  ICSH  CELL  (M  ALL  SIDES.  SINCE 
THE  DATA  BAGE  POINTS  AfC  LOCATES  ON  DC  HALF  DEGREE.  IT 
IS  BETTER  TO  PASS  VALUES  LOCATED  AT  N  ♦  1/2. 

SUHWTDC  ICSHST  IS  IN  TIC  MAIN  SIGMA  CODE 


PARAIETER  PIF*3,PJF*20 

DIMB6I0N  Z(PIF.PJF).DREE(301.91).LAM(KIF),PHI(PJF) 

REAL  *4  LAH.UMX.UIIN.UIF 

IWCLIST/DOMIN/LHAI  .LAIN.  PIMX .  PHIN.  IF  >  JF.  XLSPAC.  XPSPAC 


CALL  R4STIP 
XLSFAC-O. 

XPSPAOO. 

READ  (3.  DOMINI 

IF  (IF.BT.PIF  .OR.  JF.GT.PJF)  STOP  8GB 
PRINT  902.LNM.UMX.PMII.PmX.  IF, JF 
902  FORMAT  (1H1, '  THE  DOMAIN  RUNS  FROI  ’.F6.1 
.  .  1  KERBS  EAST  TO  '.F6.1.'  DEGREES  EAST.  AN)  ’ 

.  .FA. I.*  DEGREES  NORTH  TO  *.F6.1,’  DEORBS  NORTH.'./ 

.  .  1H0, '  THE  ICSH  IS  (IF*)  ’,13,'  BY  (JF*)  M3.'.') 

IF  (XLSPRC.B.O.)  PRINT  90S.  XLSPAC 
IF  (XP9FAC.NE.0.)  PRINT  906.  XPSPAC 
905  FORMAT  I//,*  TIE  CENTRAL  UMBITUOE  SPAC1W  IS  MPEL0.3 


906  FONMT  (//.'  TIC  (SITRAL  LATITU  K  SPACING  IS  MFEI0.3 


ROUND  t 

mam  2 

READ  (II  DREC 
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55 

CALL  ISSHGT<LAH,lF,imU«tt.IL9PAC.0> 

TOPOO 

55 

56 

CALL  lf9HBT(FHI,iF,PNIN»PHAX»XP9PAC.O) 

TQFQ6 

56 

57 

IFN1*IF  -  1 

TOPOG 

57 

58 

JFN1-JF  -  1 

TOPOO 

58 

59 

00  200  1*2.  IFNI 

TOPOG 

59 

60 

00  200  J*2„FN1 

TOPOG 

60 

61 

200  CALL  INTRPIZII.JI.LAHID.PHIIJI.OREC) 

TOPOG 

61 

62 

I 

TOPOG 

62 

63 

100  CONTINUE 

TIFOO 

63 

64 

1 

TOPOG 

64 

65 

«  FIX  U>  BOUMARIES 

TONE 

65 

66 

t 

TOPOG 

66 

67 

DO  210  tKt.JFHl 

TOP06 

67 

68 

Z<2 

TOPOG 

68 

69 

210  ZUFNl.J)*!. 

TOPOG 

69 

70 

« 

TOPOG 

70 

71 

DO  220  1*2.  IFNI 

TOPOG 

71 

72 

220  Z<I.2)*1. 

TOPOG 

72 

73 

C 

TOPOG 

73 

74 

DO  230  1*15.  IFNI 

TOPOG 

74 

75 

230  ZU.JFN1)*!. 

TONE 

75 

76 

« 

TOPOG 

76 

77 

t  00  BOMMRIES 

TOPOG 

77 

78 

* 

TOPOG 

78 

79 

DO  300  >2,JFM 

TOPOG 

79 

80 

Z(1.J)*Z(2.J) 

TOPOG 

80 

81 

300  Z(IF.J)*Z(IFN1.J) 

TOPOG 

81 

82 

DO  350  1*1.  IF 

TOPOG 

82 

83 

Z<I.1)*ZU.2) 

TOPOG 

83 

84 

350  Z(I.JF)-ZU.JFNl) 

TOPOG 

84 

85 

« 

TOPOG 

85 

86 

1  IWCE  LAW  INTO  SHALLOW  HATER 

TOPOG 

86 

87 

« 

TOPOG 

87 

88 

DO  400  iM.JF 

TOPOG 

88 

89 

00  400  1*1.  IF 

TOPOG 

89 

90 

400  IF  (Z(I.J).lE.l.)  Z(I,J>»1. 

TOPOG 

90 

91 

t 

TOPOG 

91 

92 

t  PUT  IT  UP 

TOPOG 

92 

93 

t 

TOPOG 

93 

94 

C  00  510  11*1.2 

TOPOG 

94 

95 

c  until  -  2 

TOPOG 

95 

96 

UP-2 

TOPOG 

96 

97 

WRITE  <Ui.904)LHIN.LHAX.FNIN,FNU .  IF,  JF,  XLSPAC,  XPSPAG 

TOPOG 

97 

98 

904  FORHAT  (4F7. 1,217, 1P2E10.3) 

TOPOG 

98 

99 

00  500  1*1.  IF 

TOPOG 

99 

100 

500  WRITE (LU,9Q3MZ(1,J),>1,JF> 

TOPOG 

100 

101 

903  FORHAT  (18F7.1) 

TOPOG 

101 

102 

510  CONTIWE 

TOPOG 

102 

109 

t 

TOPOG 

103 

104 

CALL  PRT 10. ’  OEPTHBA* , Z.PIF.PJF. 1, 1. .TRUE. .0) 

TOPOG 

104 

MB 

I 

TOPOO 

105 

106 

STOP 

TOPOG 

106 

107 

m 

TOPOG 

107 

108 

9M0UTIIC  INTRP  (Z.LAH.PHI.D) 

TOPOG 

108 

109 

t 

TOPOG 

109 

110 

t  THE  DATA  IS  LOCATED  AT  10*  INTERVALS.  TIC  USER  LEFT 

TOPOG 

110 

111 

I  COM  IS  AT  -10  DES  E.  AM  +30  DEG  N. 

TOPOG 

til 

112 

t 

TOPOG 

112 

E-2 
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113 

114 

115 
114 

117 

118 
114 
120 
121 
122 

123 

124 

125 
124 
127 
12B 
12? 

130 

131 

132 

133 

134 

135 
134 

137 

138 

139 

140 

141 

142 

143 

144 

145 
144 

147 

148 

149 

150 

151 

152 

153 

154 

155 
154 

157 

158 

159 


t 


« 


t 


1L  MB  X  ME  THE  INDICES  OF  TIC  NEXT  LOO  DATA  BASE  POINT.  SI 
MO  SJ  ME  TIC  FACTIONS  THAT  AC  ACTUAL  POINT  IS  ABOC  IL  AND 
JL.  THAT  IS  IF  SI-.5.  Tt4EM  AC  ACTUAL  POINT  IS  HALE  NAY  BETWEEN 
IL  AND  IL+1. 


OUB6ION  0(301.91 ) 

REALM  LAH.LOT 

DATA  LC0R/-IO./.PC0R/30./ 

XL-(LMHXOR)  t  4. 

IL-IFIHXL) 

IU-IL  ♦  1 

IF  (1L.LT.1  .OR.  1U.GT.30U  CALL  STOPP 
.('  yOA  REGION  DOES  NOT  LIE  COfLETELY  ON  DATA  SET  I') 
SI-XL  -  FLOAT!  IL) 

XMPHI-PCOR)  «  4. 

X-IFIX(XP) 

JU-JL  ♦  1 

IF  (X.LT.1  .OR.  JU.6T.91)  CAU.  STOPP 
. (’  YOUR  REGION  DOES  NOT  LIE  COFLETELY  ON  DATA  SET  *') 
SJ-XP  -  FLOAT  (X) 

SIR-1.  -  SI 
SJM.  -  SJ 

Z-  D(IL.X)  t  SIR  I  SJR 
Z-Z  ♦  D(IL.JU)  I  SI  I  SJR 
Z-Z  ♦  D(IU.X)  t  SIR  t  SJ 
Z-Z  ♦  D(IU.JU)  *  SI  I  SJ 


I 

RETURN 

END 

/  IF  E.W.O.SKIP 
/  RGET  FTOIFOOI.S/HTOPO 
/  REL  FT02F001.FT06F001.SEE 
/  FD  R02F001.LREC-124,BrSZ=2520>FOR&«.RCFMBS 
/  OXT 

•domain 

U1IIM0..  LHAX-37.. 

PNIIM2.1  PNAX-38.. 

IF-83.  JF-20. 

XLSPAOO.i  XPSPAC-O.. 


IBB 

/  PUT  FT02F001.S/ZDATA 
/  PUT  SEE.S/TI 

/SKIP  MR 


159  ACTIVE  LINE(S)  0  INACTIVE  LINE(S) 

til  96147  i  P06  DIRECTORY  FOR  tGCPL  FILE  SUCCESSFULLY  UPDATED  FOR  DECK  TOPOG 
UK  LOCATION.  ACTI9MOO 


TOPOG 

113 

TOPOG 

114 

TOPOG 

115 

TOPOG 

114 

TOPOG 

117 

TOPOG 

118 

TOPOG 

119 

TOPOG 

120 

TOPOG 

121 

TOPOG 

122 

TOPOG 

123 

AFOG 

124 

TOPOG 

125 

TOPOG 

124 

TOPOG 

127 

TOPOG 

128 

TOPOG 

129 

TOPOG 

130 

TOPOG 

131 

TOPOG 

132 

TOPOG 

133 

TOPOG 

134 

TOPOG 

135 

TOPOG 

134 

TOPOG 

137 

TOPOG 

138 

TIMS 

139 

TOPOG 

140 

TOPOG 

141 

TOPOG 

142 

TOPOG 

143 

TOPOG 

144 

TOPOG 

145 

TOPOG 

144 

TOPOG 

147 

TOPOG 

148 

TOPOG 

149 

TOPOG 

150 

TOPOG 

151 

TOPOG 

152 

TOPOG 

153 

TOPOG 

154 

TOPOG 

155 

TOPOG 

154 

TOPOG 

157 

TOPOG 

158 

TOPOG 

159 
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CREATED  ON  02/13/84  AT  11*19:04 


LAST  IFOATED  ON  02/15/84  AT  11*19*04  BY  98  VERSION  5.27 


USER  HF0ANAT10N* 

LIHIT  hiim 

F1L  LIBWRm/OBJLIB 

0IICN6I0N  XLAT  (63.63)  <  XLONB  (63.63)  <0(63.43) 

CALL  WSTOP 

wwao.  /  <4.*at«nu.» 

CALL  GRIDCU.63.1.63. XLONB. XLAT. 0) 

XLONG-XLOW  I  AD 
HAT  *XLAT  t  RD 

0-0. 

DO  10  >1.63 
00  10  1*1.63 

IF  U  XLAT<I.J).0E.32.  .AND.  XLAT<I,J).L£.3B.  )  .AW. 

.  (XLONC(I.J).GE.  10.  .AW.  XL0W(t,J).LE.37.)>  0<I,J)*1. 
10  CONTIWE 

IS-0 

15  I5=IS  ♦  1 
DO  20  >1.63 

IF  (D(ISiJ).EQ.l.l  GO  TO  25 
20  CONTINUE 
00  TO  15 
23  CONTINUE 

IE*64 

30  IE*IE  -  l 
DO  35  >1.63 

IF  (O(IE.J).ED.l.)  GO  TO  45 
35  CONTINUE 
GO  TO  30 
45  CONTINUE 

JS-0 

50  JS*J6  ♦  1 
DO  55  I*IS. IE 

IF  (D(I.J6I.EQ.t.l  GO  TO  60 
55  CONTINUE 
GO  TO  50 
60  CONTINUE 

JE*64 

65  JE-JE  -  1 
DO  70  I-IS.IE 

IF  (O(IiJE).EB.t.l  GO  TO  75 
70  CONTIWE 
00  TO  65 
75  CONTIWE 

JKJB*  JE 
DO  200  I-IS.IE 


LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATION 

LOCATIW 

LOCATION 

LOGATIOt 

LOCATION 

LOCATION 

LOCATION 


DO  90  >J8.JE 


96  VERSION  5.27  LIST1N6  OF  DECK  LOCATION 


02/15/84 


55 

IF  (D(IiJ).ED.l.)  00  TO  95 

LOCATION 

56 

90CHTIME 

LOCATION 

57 

I 

LOCATIOI 

SB 

95  CQNTME 

LOCATION 

59 

OQ  100  iMBiJE 

LOCATIOI 

60 

JL-JH-d 

LOCATION 

61 

IF  (OU.JU.ED.l.)  GOTO  105 

LOCATION 

62 

100  CQNTlMf 

LOCATION 

63 

t 

LOCATION 

64 

106  CWTINUE 

LOCATION 

65 

PRINT  902<I<JB<X 

LOCATION 

66 

902  FQRNAT  ('  AT  l>  M2.*  J  RUNS  FRDH  M2.'  TO  M2) 

LOCATION 

67 

200  CONTINUE 

LOCATION 

68 

1 

LOCATION 

69 

1  LEAVE  AN  EXTRA  (ME  ALL  TT€  MAY  ARMS  FOR  T>€  INTERPOLATION 

LOCATION 

70 

t 

LOCATION 

71 

IS-IS  -  1 

LOCATION 

72 

IE-IE  ♦  1 

LOCATION 

73 

JS-JS  -  1 

LOCATION 

74 

JE«X  ♦  1 

LOCATIOI 

75 

t 

LOCATION 

76 

PRINT  901.IS.IE.JS.JE 

LOCATION 

77 

901  FORNAT  <*  I9»M2.'  !E*M2.'  J9*M2,'  JE«M2> 

LOCATION 

78 

EM) 

LOCATION 

79 

/  OXT 

LOCATION 

80 

/  CON  ANSWER: 

LOCATION 

80  ACTIVE  LIIC(S)  0  INACTIVE  LIIE(S> 

96167  :  PCS  DIRECTORY  FOR  NEMGPl  FILE  SUCCESSFULLY  UPDATED  FOR  DECK  LOCATION 

**  INITIAL.  ACTIOHAOB 


»  VERSION  S.27  LISTING  OF  DECK  INITIAL 


02/15/84 


CREATED  ON  03/19/M  AT  11*1*08  LAST  MATED  ON  02/19/M  AT  11*  19*09  IV  SHB  VERSION  9.27 

LANBUMEi  USB)  IMRHATION* 


1 

/  SET  08XIB.PT/aULIB 

initi 

2 

/  UK 

INITI 

3 

LIBRARY  OULU 

INITI 

4 

INCLUDE  TAFER3 

INITI 

5 

/  TAPEAEAD  9EE0PT-O,  TAPEPATHM/CL I NAT.  CATPATW-S/  INI TUTTA,  i 

INITI 

& 

IWIIC«2000.FILESIZE«20.L0AW»* 

INITI 

7 

me 

INITI 

8 

ISMS.  IE-90.  JS-31.  JEMO. 

INITI 

9 

M>20.  SIONAMT. 

INITI 

10 

IBB 

INITI 

11 

T  SEAA 

INITI 

12 

T0400  A 

INITI 

13 

T0600  A 

INITI 

14 

TOBOO  A 

INITI 

19 

T1000  A 

INITI 

16 

T1S00  A 

INITI 

17 

T2000  A 

INITI 

18 

T3000  A 

INITI 

19 

T4000  A 

INITI 

20 

T9000  A 

INITI 

21 

SOOOO  A 

INITI 

22 

90090  A 

INITI 

23 

90100  A 

INITI 

24 

S0200  A 

INITI 

29 

90600  A 

INITI 

26 

S1000  A 

INITI 

Z7 

92000  A 

INITI 

28 

S3000  A 

INITI 

29 

94000  A 

INITI 

30 

99000  A 

INITI 

30  ACTIVE  LIIEIS)  0  INACTIVE  L1NE(S) 


*11  9*8167  *  PD6  DIRECTORY  HR  IBCPL  FILE  9UCCES9ULY  UPDATED  FOR  DECK  INITIAL 


SHE  VERSION  9.27  LISTING  OF  DECK  FORCING 


02/19/84 


CRBUEB  ON  02/19/84  AT  lit  1*06  LAST  UPDATED  ON  02/19/84  AT  ill  1*06  BY  SHB  VERSION  9.27 

LANBUMEi  USER  MFORHATION 


1 

/  ABET  OBJLID.pt/OBXIB 

FORCING 

2 

/  UK 

F0CIN6 

3 

LIBRARY  QBJLIB 

FORCING 

4 

INCLUDE  TAPER 

FORCING 

9 

/  TAPEREAD  TAPEPAT)69(/SAIATTC.CATPAT>#S/FCRCI*. ! 

FORCING 

6 

RUNTIfE-4000,FILESIZE«20,UDW«t 

FORCING 

7 

MNE 

F0RC1IC 

8 

YEAR-77,  NKDM.  DAY-7,  NOUN), 

FORCING 

9 

19-45.  IE-50,  JS-31,  JE-40, 

FORCIIG 

10 

NT-96.  DT-6., 

FORCING 

11 

SIGHAD-T,  A27A»T, 

FORCING 

12 

IBS 

FORCING 

12  ACTIVE  LINE(S)  0  INACTIVE  LlfC(S) 


Ut  9HB167  i  PUS  DIRECTORY  FOR  KUSPL  FILE  SUCCESSFULLY  UPDATED  FOR  DECK  FORCING 


